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We study the weakly non- linear development of shear-driven gravity waves, following 
the physical mechanism first proposed by Miles and furthermore investigate the mixing 
properties of the finite amplitude solutions. Calculations to date have been restricted 
to the linear theory, which predicts that gravity waves are amplified by an influx of en- 
ergy through the critical layer, where the velocity of the wind equals the wave phase 
velocity. Because of the presence of a critical layer, ordinary weakly non-linear meth- 
ods fail; in this paper, we use a rescaling at the critical layer and matched asymptotics 
to derive an amplitude equation for the most unstable wave, under the simplifying as- 
sumption that the physical domain is periodic. These amplitude equations are solved 
numerically, in their quasi-steady limit, for the cases of small density ratio (applicable 
to oceanography), and for arbitrary density ratio but strong stratification (for more gen- 
eral physical/ astrophysical situations). As is found in other analysis for critical layers in 
inviscid parallel flow, we find that the initial exponential increase of the amplitude A 
transitions to an algebraic growth rate proportional to the viscosity, A ~ vt^^^. How- 
ever, for the air over water case (provided the maximum wind velocity is in the range 
of 0.2 ms~^ ~ 1 ms~^), our results from the weakly non-linear analysis for the single 
mode show that the transition from exponential to algebraic growth rate occurs when 
the amplitude of the wave is as small as h ^ lO^^A; hence, it may be difficult to observe 
the linear regime for this case in numerical simulations. (Whether this result transfers to 
cases in which the physical domain is not periodic, so that the marginally stable modes 
form a continuum, is not as yet known.) We also find that the weakly non- linear flow 
allows for super-diffusive particle transport with an exponent ~ 3/2, consistent with 
Venkataramani's results. 
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1. Introduction 

The generation of surface waves by winds has been a problem under study for weU over 
a century. The simple model of Kelvin-Helmholtz (from now on, KH) instability provided 
higher bounds on the maximum wind velocity for the instability to occur than what was 



experimentally observed (Chandrasekhar (1961), §101); and a solution was not available 



until Miles (1957, 1959a b 1962 1967) proposed a model in which gravity ocean waves are 
amplified through a 'resonant' interaction with the wind above the ocean surface: Waves 
that travel with speed c = \/g/k are amplified through an influx of energy at the height 
Uc in the wind where the wind speed Uijj) matches the wave speed, U{y = yc) = c. Miles 
explicitly assumed that the wind profile above the sea surface was of the form log(z/z*) 
(where z* is a parameter), composed of a mean velocity U and a turbulent component 
u. Assuming the mean wind profile to be given from turbulent boundary layer theory 
and the turbulent component to be small, he calculated the influx of energy to the wave. 
This analysis involved the treatment of the critical layer (the point where the velocity 
of the wind was equal to the velocity of the gravity wave), where the solution of the 
linear eigenfunction problem becomes singular (the x-component of the perturbation 
velocity behaves as u ~ log(y — yc))- The singularity is removed either because there is 
an imaginary component of c (i.e., unsteady critical layer) or because viscosity becomes 
important in this thin layer (i.e., viscous critical layer). In both cases, the final result 
is that there is an '— itt phase change' across the critical layer, e.g., the perturbation 
wave above the critical layer is not in phase with the wave below. A direct result of this 
phase change is that the gravity wave will become unstable because a component of the 
pressure perturbation will be in phase with the slope of the wave (unlike the KH case, 
in which the pressure perturbation wave is always in phase with the wave). Using the 
results from linear theory and for small density ratios. Miles predicted the energy flux 
from the wind to the gravity wave. 

Our interest in this problem is motivated by an astrophysical puzzle, namely, the 
mixing of carbon/oxygen (C/0) at the surface of a white dwarf with accreted material 
(mostly hydrogen and helium, He/H) overlying the stellar surface (whose presumed origin 
is from an accretion disk surrounding the compact star). For a variety of reasons, it is 
thought that the accreted envelope is in differential rotation with respect to the stellar 



surface, so that a 'wind' is expected at this surface (Rosner, Alexakis, Young fc Truran 



Hillebrandt 2001 ). In this stellar case, one is forced to generalize the earlier results to 
arbitrary density ratios (between the 'atmosphere' and the 'surface material'); and we 
have already done so for the linear problem ( Alexakis, Young &: Rosner 2002| ), deriving 



bounds on the instability in the parameter space, and estimating the growth rates of 
the unstable modes. While useful as a starting point for further analysis, linear theory 
gives no information about the ultimate fate of the system, which is largely governed by 
nonlinear processes. In this paper, we uncover the non-linear effects on the evolution of 
wind-driven surface waves by examining their finite-amplitude evolution. To explore the 
applications to the relevant astrophysical and geophysical systems, we also investigate 
the mixing property of the weakly non-linear flow using particle method, and results 
show enhanced mixing due to the coupling of the weakly nonlinear wind with the surface 
wave. 

We note here that the weakly non-linear theory for the KH instability has already 



been derived by Drazin (197C). That analysis cannot be applied straightforwardly to our 
problem because of the presence of critical layers: Due to the singularity that appear in 
the linear theory at the point where the phase speed of a surface wave matches the wind 
speed, higher order terms in the expansion become more singular and the expansion must 
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ultimately fail. The fundamental reason for this behavior is that the flow becomes non- 
linear first inside the critical layer even though the rest of the flow can still be considered 
as operating in the linear regime. For this reason, a more refined treatment of the critical 
layer is required. The necessary analytical 'machinery' fortunately already exists: Thus, 



Benney & Bergon (1969), and later Benney & Maslowe (1975), observed that for small 



but finite amplitudes the phase change at the critical layer does not necessarily have to 
be —iTT; instead, they found a solution with the property that if the nonlinear terms are 
taken into account, the phase change is zero. Later Haberman (1972) showed numeri- 
cally that there is a smooth increase of the phase change, from —in to 0, and introduced 



the function <& that gives th e phase change as a function of the amplitude; phurilov fc 
Shukman (1987 , 1988| , 1996 ) then developed a weakly nonlinear theory based on $ and 



other similar functions defined for the appropriate critical layer problem. A fundamental 
assumption in all this work is that the viscosity is dominant in the critical layer; this 
leads to the derivation of an ordinary differential equation for the wave amplitude. The 
full equations of the weakly nonlinear problem without making the previous assumption 
have then been solved numerically for various cases (Goldstein fc Hultgen 198S; Warn fc 



Warn 1978| ; [Balmforth fc Piccolo 20011 ). 

In our case the treatment - although closely related to the earlier work - is nevertheless 
different in two respects. The first one is that the previous cases possess stability bound- 
aries for modes for which the critical layer is formed at the inflection point dyll/dy^ = 0; 
one then examines the close-to-marginally-stable state unstable modes. This is however 
not our case since such an inflection point does not exist (unless one considers C/(oo) as 
such a point). Instead, we will be considering two special cases: the case of small density 
ratio, and the case of very strong stratiflcation^; in both cases, the linear growth rate is 
small. This procedure will lead to a slightly different scaling. The second difference lies 
in the fact that we have an interface. Because of this, our solvability condition will not 
be expressed in terms of integrals but rather in terms of appropriate vector products of 
the values of the perturbation stream function at the interface. 

This paper is structured as follows. First we formulate our problem and describe the 
non-dimensionalized form of the basic problem equati ons. In §3 we very briefly review 
(and slightly revise) the linear theory developed earlier ( Alexakis, Young fc Rosner 2002 ). 
The amplitude equations are derived in §4.1 and §4.2 for the small density ratio and 
the strongly stratified cases, respectively (we define these terms more precisely in what 
follows). We summarize the conservation laws from the amplitude equations in 
and we discuss the implications of the long-time behavior of surface waves in 



5.1 



§5.2 



In 

we summarize results from numerically simulating the amplitude equations. In §6 
we investigate how the coupling between the wind and the surface wave affects mixing 
properties of the finite amplitude solution. A detailed examination of the assumptions 
made in the analysis is provided in §7, where we also draw conclusions from our work. 



2. Formulation 

We consider a two-layer system with constant fluid density pi and p2 {pi < P2), 
respectively, in the upper and lower layers. The interface between the two layers is given 
by y = h{x,t), where x is the horizontal and y is the vertical coordinate. In the upper 
layer we assume a wind parallel to the originally flat interface y = h(x,t = 0) = 0. The 

f This case is of particular interest to our astrophysical problem, since the instability takes 
place on the surface of a white dwarf star, i.e., a star of solar mass, but comparable in size to 
the Earth. 
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wind has a shear velocity profile U = [U{y), 0], where U{y) varies only with y. The lower 
layer is initially at rest. Our aim is to study the dynamics and the weakly nonlinear 
development of a small sinusoidal perturbations of the horizontal interface.^ 

The fluid is assumed incompressible in both layers; we therefore work with the stream 
function 'J, which is connected to the velocity by the relation {u, v) — — The 

stream function can be separated into its mean and its perturbation components, 

*totai= Cudy' + -^. (2.1) 

JO 

The fluid within each layer is described by the Navier-Stokes equations (in terms of the 
stream function ^E"): 



VH^t + UVH,, - U^yy^^^ = ^,=.V2*,j, - *,aV2*^, + vV^VH, (2.2) 

where V is the two-dimensional Laplacian and v is the viscosity (which we will consider 
to be small, and therefore negligible except for a narrow region within the critical layer). 
Here we have used the standard notational device of comma-prefaced subscripts to denote 
partial derivatives, e.g., 

f,x = df/dx. 

The boundary conditions at the interface are the continuity of the perpendicular com- 
ponent of the velocity to the interface 



+ (C/± + *±)/i,, + 'J'^ - , (2.3) 

and the continuity of pressure 

,yy) 

= gh.,x{P2- Pi), (2.4) 
where h{x) is the elevation of the interface and all quantities above are evaluated at 
y = h(x). The ± indices indicate values above and below the interface, and A[ ] denotes 
the difference across y — h (e.g., A[/(?/)] = f{h'^) — f{h'~) ). 

We non-dimensionalize lengths by the characteristic length I of the wind, and velocities 
by the asymptotic value of the wind at y ^ +oo, Umax- An important parameter that 
emerges from the scaling v& G = gljU^ , which is a measure of the ratio of potential energy 
to kinetic energy, or alternatively, a measure of the strength of the stratification. Other 
dimensionless parameters are the Reynolds number. Re = Uljv ^ 1, and the ratio of 
densities, r = pi/ P2 < 1- In terms of the non-dimensional parameters, we focus on cases 

X As an aside, we note that the presence of a sharp interface boundary between the two fluids 
in our system is really only a device to simplify the analysis, but in no way restricts our results. 
That is, in a more general case there will be no sharp interface, but the density will change 
smoothly from p2 to pi in a layer of width 5\. We would then also expect U(y) to be a smooth 
function in y, so that any flow discontinuity at the interface would be replaced by a smooth 
variation in U within a thin viscous boundary layer of width 52- One would then expect to see 
differences in behavior only for those modes with horizontal wavenumber k large enough so that 
the critical layer lies inside these layers 5i,52- Such modes, however, are KH-modes which are 
not under study here; they will obey a Richardson-type stability criterion, 1/4 < g{Ap/ p)52/S 



( Chandrasekhar 1961). We therefore expect our assumption of a sharp interface to be reasonable 
for horizontal wave-numbers > max{5i, ^2}. Since we are primarily interested in long wave- 
length perturbations, we conclude that, for our purposes, we have not disregarded any important 
physics. 
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where Re ^ 1 in each layer, and especially the two distinguished limits r <^ 1 (air on 
water, for example) and arbitrary 'large' density ratio r with G 3> 1 (accretion on white 
dwarfs in astrophysics). 



3. Linear Theory 



We linearize equations (2.2) - (2.4), and write the stream function perturbation ^> and 



surface elevation ft as a traveling wave parallel to the wind (right traveling wave in our 
setup) 



vl>±(y,x,t) = 0±(y)e^^(^-^*)+c.c., 
h = /ie*^(--<^*) + c.c. , 

where c.c. denotes the complex conjugate, and K — kl is the non-dimensional (horizontal) 
wavenumber and C is the non-dimensional phase velocity. This leads us to the well- 



studied Rayleigh equation (Drazin & Reid 1981) 



^,yy 



U 



yy 



= 0, 



(3.1) 



u ~c 

where the inviscid limit Re ^ oo is taken. The boundary conditions become 

iU^ -C)h + 4)^ (3.2) 

A[p,(({/±-C)0,,-C/±0±)] -G/i(p2-Pi)=0, (3.3) 

Eliminating 0" and h from the last two equations, and dropping the + index for conve- 
nience, we obtain 



KC^ -r 



where we have used 



{U -Cfc^,y^[U -C)U,y -G(l-r) = 0, 



1 as a normalization condition. 



(3.4) 



To further simplify we assume that U{Q) — and omit the KH- modes (of which the 
nonlinear evolution is not under study here, see Alexakis, Young fc Rosner (2002D for 
more details). The linear growth rate (KlrcL{C) ) has been numerically calculated and 
summarized in Alexakis, Young fc Rosner (2002 ) for cases of interests. Here we briefly 



summarize some relevant results on stability boundaries for the wind profile 



U = 1 



-y . (3.5) 

As shown in Alexakis, Young & Rosner (2002), this wind profile allows for an analytic 
expression for the stability boundaries in the stability diagram. One can show that some 
general features of the stability boundaries summarized here for wind profile in equation 
( p.5[ ) also hold for other bounded wind profiles. The stability bound comes from the 
modes that have phase velocity G — 1, in which case the solution of equation (3.1) 
becomes — e***^ with k = Vl + ■ Applying the boundary conditions, one obtains 
the criterion on wavenumber K for the unstable modes: 



K > K„ 



G{l-r) + r- r^(G(l-r)+r)2 + (l-r2) 



1 



(3.6) 



The unstable modes can be further restricted if we assume the presence of surface tension 
or magnetic fields. This will lead to an additional term TK"^ in equation (3.4), where 
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T = o'/(/32C^maa;0 ^"^d CT is the surface tensiori|f|. Then the unstable modes he in the region 
Kmin < K < Kmax, where Krnin and Kmax are given by the positive real solutions of 

K + rVl + K^- (G(l -r)+r + T{K)K'^) = . (3.7) 
Furthermore, there is a minimum value of T below which the above equation has no 
positive real solutions, and therefore no unstable modes exist. The physics behind these 
bounds is simple: In order for a mode to become unstable for the above wind profile, the 



phase velocity of the wave must lie in the range < c < Umax ( Alexakis, Young fc Rosner 



2002 ). With the inclusion of surface tension, the phase velocity is not monotonic with 
K , but for large enough K , the phase velocity increases in an unbounded fashion. This 
leaves only a finite region in K space with phase velocity smaller than Umax] moreover, 
if the surface tension is large enough, the minimum phase velocity is larger than Umax, 
and therefore no unstable mode exists. 

The analysis is simplified if we assume small density ratio, and expand all quantities 
in r: 

C = Co + rCi + . . . 
(/) 00 + f 01 + • • ■ 

In that case the linear theory, to zero order, gives a gravity wave with phase velocity 
Co = yjGjK. At the next order, one obtains 

2CoCi = [C^^o.a + V^y\ - G , (3.8) 



and therefore 



where Sn is such that 



Im{Ci} = icim{0o, J (3.9) 



00,,,- ^' + 77%r -^0 = 0. (3.10) 
(7 — Oo_ 

An —ITT phase change at the critical level is assumed. Due to this phase change, 0o.j/ is 
complex at the interface and thus C\ has an imaginary component, which is responsible 
for the instability of the traveling waves. 

For the arbitrary r case, the previous expansion does not hold, and one needs to solve 
the full set of equations (3.1)-(3.4) as in Alexakis, Young fc Rosner (2002| ). There it was 



shown that the growth rate exhibits an exponential dependence on the parameter G. In 
Appendix A, we carry out an asymptotic analysis for large G, and derive this exponential 
dependence. More specifically, the growth rate is found to be 



K Im{C} 



1 



4(1 - r) AtG 



3/2 



(KjAtG) 



^K/{AtG) 



2 AtG 



(3.11) 



where At is the Atwood number. As is shown in the Appendix, the stream function 
below the critical layer is composed of an exponentially increasing and an exponentially 
decreasing component. Since the boundary condition (equation ( |3.2| )) must be satisfied, 
the exponentially large component must be in phase with h, which leaves us with the 
exponentially small, out of phase, component to drive the wave unstable. A detailed 
calculation then leads to the result given above. 



t We note that a magnetic field whose direction is al igned with the interf ace will have the 
same effect as surface tension with (j{K) = ( Chandrasekhar 1961). 
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4. Weakly nonlinear theory 

We are now ready to embark upon the weakly non-linear theory. Formally this is 
done by assuming that, for some parameter ranges of interests, our physical system lies 
close to a marginally stable state so that an asymptotic expansion is allowed near the 
center manifold. For the problem at hand, however, and in the absence of surface tension 
(T = 0), marginally stable states are possible only when r = Oorl/G' = (we note 
that K is not a control parameter); the first one expresses the unphysical situation that 
there is no upper fluid, and the latter corresponds to a situation where there is no wind 
in the upper fluid {U ^ 1/VG)- Complication arises as we deviate from these neutrally 
stable states. In ordinary dissipative systems, only a small number of modes near the 
center manifold become unstable and need to be considered. In our case though, once 
the density ratio r or the parameter G is finite, an infinite number of modes become 
unstable if surface tension T — 0. Ideally the interaction of all these modes needs to be 
taken into account. Practically this difficulty is removed by the combined effect of surface 
tension T (or, equivalently, the presence of a magnetic field) and weak viscous damping. 
In the presence of surface tension the stability boundary in the (r, G, T) space is given 



by the condition for positive solutions of eq.(|3.7|) for r = 0. Surface tension reduces the 
number of unstable modes by neutralizing modes of wave numbers above some cut-off 
value. Further more weak viscosity will damp out the neutrally stable modes of high 
wave numbers, rendering them asymptotically stable. Using a periodic domain we can 
always fix the period so that only one mode becomes unstable. We can then derive the 
amplitude equation for this single mode as usual. 

We find that these amplitude equations are almost the same as those obtained from 
asymptotic expansion around the most unstable mode without surface tension; the only 
difference is a coefficient dependence on the surface tension. In the following derivation, 
for simplicity, we chose not to include surface tension, and assume the presence of surface 
tension only implicitly. At the end of the analysis the surface tension can always be 
recovered by replacing G with G + TK^ . 

With the previously stated assumptions, and following the basic strategy of the weakly 
nonlinear theory, we develop below an asymptotic expansion based on the small amplitude 
of the perturbation, and introduce a long time scale of the same order as the nonlinear 
terms inside the critical layer. We consider two different cases: The first case assumes that 
the density ratio is small (r ^ 1); the second case allows for arbitrary r but assumes that 
G ^ 1. The samll density ratio (r) in the first case is the parameter to which we scale 
our amplitudes. The large G and arbitrary r case, however, is a bit more complicated; 
for this case, the growth rate according to linear theory is proportional to exp(— 2G), 
and this will be our scaling parameter. 

4.1. Small density ratio case 

We start with the small r case. The following analysis applies to for a monotonicly 
increasing but otherwise arbitrary wind-profile. The outline of the derivation is as follows: 
we write the Euler equations, and the boundary conditions, in the reference frame of the 
moving wave. We introduce the scaling discussed below, and solve up to the second order; 
this requires a more detailed treatment at the critical layer that determines the phase 
change across the layer. Using matched asymptotics, we match the inner solution with 
the Frobenius solutions (j)a, 4>b of the outer flow. This introduces five different amplitudes, 
, , and A2 , for each Frobenius solution above and below the critical layer (which 
will be assumed to be given numerically); the fifth amplitude A2 is just the surface wave 
amplitude. The first four amplitudes are connected through the boundary conditions at 
infinity and the matching at the critical layer; the fifth is connected with the rest through 
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the boundary conditions at the interface. When we apply the boundary conditions, the 
first order gives the free gravity wave (vacuum over water); at second order, the solvability 
condition results in the amplitude equation. 

We begin by introducing the scaling. Let e = pi/p2 1 and dt = cOt — Cdx- Then 
the equation for the stream function in the reference frame of the wave will be: 

eVH^T + {U- C)VH^, - U^yy^^, = *,xV2*,j, - 4-,^ V^^-,, + Iv^V^* , (4.1) 
with boundary conditions at y = 

eh^T — Ch^x + = non-linear terms (4.2) 



and 



Gh^xi^ — ^) = non-linear terms , 

(4.3) 

where the ± index means above or below the interface. First we focus on the upper fluid 
and the outer solution (outside the critical layer). 

4.1.1. Outer solution in upper fluid 
We expand the stream function as 

^ = e^^o + e^*i + • • • (4.4) 

To zeroth order we have 

Focusing on the most unstable mode and writing '^q = A{T)(l){y)e^^^ -\- c.c, we obtain: 



2 I '-'.yy 



= 0. 



u -c 

Near the critical layer y ^ Vc the solution can be written as 

The ± signs correspond to above or below the critical layer, and the two Frobenius 
solutions are 



(4.5) 



^« = (^l + ^(y-yc)ln|2/-yc| + -^ , 

(^(2/- 2/c) + ...) , (4.6) 

where U'^ and /7"are dU / dy\y=y^ and d'^U/dy'^\y^y^, respectively. In general, and 
are functions of T, so that it is more convenient to write the solution as 

*o= [^i±(r)'/'a(y) + Si±(r)</.b(y)]e'^" + c.c. (4.7) 
At the next order, we have 



Wind-driven gravity waves 
Again expanding ^'i in terms of Frobenius solutions, we write ^'i as 



iKx 



where Ai± and Bi± are the ampUtudes of ^'i. In order to match with the inner solution, 
we require the behavior of ^' as (?/ — yc) eY . Upon expansion in Y we obtain 



^ = e *o + e *i + ••• 



U" U" 
A^±-^Y^dTA^± " 



c/;2 



U" U" U" 

Bi±^Y + Ai±-i^y In |r| + arAii-^ In |r| 



e*^^^ + c.c. 



4.1.2. Inner solution 

To capture the dynamics inside the critical layer, we have to use the scaling ^ 



e^^iY), y — yc ^ eY and 1/R ^ e^v. From equation (4.1) we then obtain 

^,TYY + U^Y^^^YY + *,y*.Yyx - ^^x'^^YYY - ,YYYY 

"1 



-U'^Y''^,YYx - U'^^,, 



In order to match with the outer solution wc expand as 

* = §0 + eln(e)#i + e§2 + ln(e)^3 + e^§4 + •- • 
To first order, we then have 

^0,TYY + U'cY^Q^xYY + ^ox'^oYYx - ^0,x^O,YYY - V^O^YYYY = . 

Matching with the outer solution we obtain 



*o = Ai+e 



iKx 



A*e 



-iKx 



and therefore 



(4.8) 



(4.9) 



Ai+ = Ai_ = Ai . 

To second order (eln(e)), we have 

^i,TYY + U'jifi^^YY - ^o^x^i.yyy - i^^i.yyyy = . 
Matching with the outer solution wc obtain 

TTll 1 

AKx 



(4.10) 



U" U" 



+ c.c. 



To third order (e^), we have 

*2,yyT + U'cY^2,YYx - *o,2:*2,yyy - i^*2,yyyy = U'^^qx ■ 
Denoting Z = 5*2, yy, which is the vorticity inside the critical layer, we obtain 

Z^T + KYZ^, - ^0,xZ,Y - vZyY = t/"«'0,x • 

To match with the outer solution, we require the boundary conditions of Z 



(4.11) 



(4.12) 



(4.13) 
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lim Z 



. u'j 1 u'j 1 - 



e^-^^ + c.c. 



lim y = e 

^±00 



.Kx 



U" U" 1 U" 

Ul. ^ ' U'Y U' 



Integrating Z along Y , we obtain 



r+oo 




+ 00 


1 ZdY = 






J —oo 




— oo 



u. 



where in the above integral we assumed the limiting procedure lime^o /jti/e ^ dY . Let 



J = 7rTTJj e"'""" / ^^^^^ = ^1+ - ^1- 

iTT U'J^ J-n/K J -co 



(4.14) 



This leads to the important result: 



(4.15) 



The last equation ( 4.15 ) implies that the phase change across the critical layer depends 
on the more detailed treatment of the vorticity dynamics inside the critical layer. In the 
linear case J = +iTTAi, but as the amplitude grows and the non-linear term ^q^xZ.y in 
equation (4.12) becomes important, the phase change is going to decrease. 



4.1.3. Lower fluid 

Using the same scaling as for the outer solution of the upper fluid, we have 



To first order, we then have 



A2(T)0(7e'^"+^^^ + c., 



At next order. 



= A^(T)<^^e'^"+^^ + c.c. 

At this point we are ready to apply the boundary conditions at the interface, and obtain 
the amplitude equation. 

4.1.4. Boundary conditions and the amplitude equation 
First, we expand the amplitude of the perturbed interface 

h = e^ho + e^hi + ... 

To first order, we obtain 

ho = A2hoe^^'' + c.c. 



Then using equations (|4.2| ) and (^^), we obtain 



4>Q - Cho = , 



(4.16) 
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+ G/Kho . 



(4.17) 



Let 



and M = 



-C 1 
G/K -C 



Then the previous set of equations ( 4.16 )-( [t.l7l ) can be written as MVo = 0. For a 
non-trivial solution we must have det[M] = 0; from this condition we obtain 



and 



From equation (E^) we also obtain 



or 



C=^G/K, 
4>Q = Cho . 



(4.18) 
(4.19) 



To the next order we have 
and 



- ChoA2 + (Ai0a|o + Bi^cjjblo) = . 
hi = A2(T)hie'^'' + c.c. 



^i,x\o - C'hi.x = -hQ,T , 



or 



(4.20) 



01 ~ Chi = ~—A2,Tho , 
iK 

-C4>i + {g/k)hi - —La2^t4>o + MiG/K) - U,yC{l/k))ho - lc{Ai(ba,y + Bi_4>b,y) 



Set Vi = [hi,^i]'^ and 



Wi = 







(G/K) ~ UyC/K 



, W2 = 



-C/Kcpa.,y 





-C/K(j3h^y 



then the previous set of equations can be written as 

MVi = — -A2 tVo + A2W1 + A1W2 + Bi^W^ . 
iK 

Let V"^ = [C,l]. V'^ is chosen such that V^M = 0. Then multiplying the previous 
equation with V"^ yields 



1 

ik 



where 



■^A2t{V^Vo) = A2{V^Wi) + Ai{V^W2) + Bi^iV^Ws) 



{V^Vo) = Cho + <Po = 2Cho =- 2C , 
(V^Wi) = {G/K)-U,yC/K, 



(4.21) 



12 Alexakis, Young and Rosner 

{V^W2) = -C/K^a.y , 

where with no loss of generahty we have set hg ~ 1; here 
of these quantities at the interface. 

Equation (4.21) is the amplitude equation. What is left to do is to recall the properties 



refers to the values 



of A2 , Af , B^in order to cast the amplitude equation into a more convenient form. From 
the boundary condition at infinity, we will obtain a condition of the form = /LtAi+, 
e.g., only one linear combination of the two Frobe nius solutions will give a solution that 
dies at infinity. From equati ons (4.10 ) and ( 4.15 ) we have Ai+ = Ai- = Ai, Bi- = 
Bi+ — J, and from equation ( 4.20 ) we obtain —CA2 + {Ai(j)a + -Bi_(/>b) = 0. Using these 
results, we can write 



Ai 



Bi- 



C 



-A2 



-A2 



-J, 



<t>a 



We can therefore write the amplitude equation as 



-J. 



(4.22) 
(4.23) 



A 



2T 



(U,y - G/C) + C 



4'a,y + ^J■4'b,y 



A2 = i 



T>b(Pa,y 



Va9b,y 



J 



(4.24) 



4.2. Strong gravitation and arbitrary density ratio: G 3> 1 and < r < 1 

We are now going to focus on cases for which G ^ 1, which are associated with small 
growth rate. As before, we use the time scaling dt — edr — Cdx, where e is of the order 
of the linear growth rate (to be defined below). The governing equation is then equation 
(4.1), with the boundary condition now given by 



eh^T - Ch,.^ + = , 
eh.T - Ch,a: + ^-^^ = 0, 



(4.25) 
(4.26) 



^ ,Ty ^^,xy 



Gh,xil - r) = 0. 



We expand the surface elevation and stream function as: 

* = e^fo + + ■ • ■ 

h = e^/io + e^hi + ... 
To zeroth order we then have again 

expanding in normal modes, ^0 = ^20(y)e*^^ + c.c, we have 



0; 



U-C 



= 0, 



(4.27) 

(4.28) 
(4.29) 



where we have again focused on the most unstable mode. Since we are dealing with the 
large G case, we know from equation ( |3.6[ ) that the wave-numbers of the unstable modes 
must satisfy (in the limit of large G) 
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K > 



1 



-G = AG > 1 



(4.30) 



1 + r 

(where At is the Atwood number); we can then use the resuhs from the WKBJ expansion 
discussed in Appendix A. These results show that the value of \E'o and ^'o,y at the interface 
is 



Ai(sin(/i)e~^^ +cos(Ji)e+^=^) - -Jsin(/i) 

TT 



*0, 



v\o 



Ai(sin(/i)e~-^^ - cos(/i)e+^2) _|_ i Jsin(/i)e+^=^ 



e*^(^-^*) + c.c. (4.31) 
e'^^'^-C') + c.c. , (4.32) 



where Ii ~ 1/K and I2 ^ K > AtG ^ 1 are defined in the Appen dix. U nlike the linear 
case, we assumed that the phase change is J, defined as in equation (4.14) in the previous 
section (§4.1.2). We avoid repeating the inner scale calculations since they are identical 
to that for the small density ratio case. The slow time scale is now defined by the value 
of C'i, which is exponentially small and is given by expression (3.11). We will therefore 
define e = e^^^^. Equations (4.31) and (4.32) can then be written as 



and 



-K 



Ai cos(/i) Jsin(/i) + eAi sin(/i) 



Ai cos(/i) Jsin(/i) — eAi sin(Gi) 



c.c. . 



(4.33) 
(4.34) 



where we have rescaled \l/o so that it is of order unity at the interface. 

Writing the zeroth order stream function ^'q below the interface, and the surface 
elevation ho as 



7 A iKx I 

fto = A2e + 
the first order boundary conditions give us 

M • = 

where 



c.c. 



(4.35) 
(4.36) 

(4.37) 



-c 1 1 r A2 ■ 

M = -C 1 , = $ , (4.38) 

-G(l - r) rCK CK \ [ _ 

and $ = (2A2Cos(/i) — ijsin(/i)) is the amplitude of the stream function at the in- 
terface. Here we have discarded the Uy term since it is of order 1/K. For a non-trivial 
solution, we must have det(M) = 0. This leads us to the relations 



and 



G = ^/GAt/K 

A2G= $, 

A2C = A3 . 



(4.39) 

(4.40) 
(4.41) 
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At the next order, we have 

and by letting vj/i = Aie*^^0i(t/), we have 



u ~c 



2 ^-.VV 



(4.42) 



Since we are still in the large K limit, we can still use the WKBJ approximation for the 
above non-homogeneous equation. One needs to notice that the non-homogeneous term 
is going to be of order \/K^ everywhere except close to the critical layer, and that special 
care is needed to be taken there. Fortunately, the exact value of 0i at the interface is not 
needed — it is sufficient to know that the solution of equation ( 4.42| ) below and away 
from the critical layer is 

^^^le+^/o"^"''^'e"^« + il,le-^/o"""^'e+^^ (4.43) 
w w 

where w = {K^ — U^yy/{U — C))^/^, and where the values of Aia,Ait, depend on Ait, 
and can be obtained by matching with the inner solution at the critical layer. The second 
term in equation ( 4.43| ) is exponentially small (order e), so that at the interface we have 
(j)i^y = —K(j)i; this suffices for the remainder of our analysis. The stream function of the 
lower fluid 4"^ and the elevation hi at this order are 



^^^rKx+Ky _^ ^^^^ 



hi = A2e'^'' +C.C 



The boundary conditions at the interface are 

1 



M ■ Vi 



iK 



drWi + W2 ; 



(4.44) 



(where M was defined in equation 4. 38) 





' A2 - 




A2 




1 


Vi = 


(piAi 


,Wi = 


A2 


= A2 


1 




. As _ 




-X(r$ + A3) 




-{r + l)KC 



and 



(4.45) 



-sin(/i)Ai 

W2= . (4.46) 

rKCsm{Ii)Ai 

As before, let V'^ = [-rCK, -CK, 1] where V'^ is such that V'^M = 0. Muhiplying 
equation (4.44) with V'^ , we obtain 

V^Wi = -2(1 + r)CKA2, V^W2 = 2rCKA2 sin(/i) . 
Hence, we arrive at the amplitude equation for large G 



1 T 
iK 



or 



A2.T = ~iK——svcL{Ii)Ai , 



(4.47) 
(4.48) 
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with 



CA2 = 2Ai cos(/i) - J- sin(/i) . 



(4.49) 



The ampHtude evolution can now be determined from equations ( 4.48 )-( 4.49 ) 



5. Results 

5.1. Preliminaries 

For both cases we have considered, the ampUtude equations to be solved in order to 
determine the nonlinear evolution of the unstable mode can be written in to the following 
more general form: 



A2.T + iCiA2 = J 



(5.1) 



A2 ^ViAi+V2J 



(5.2) 



2tt U'J 



-iKx 



+00 



Z dYdx = 



iKx 



z) 



Z^T + U'JZ,, - ^o,xZ,Y - vZyy = C/;§o,: 



(5.3) 



(5.4) 



with 



and 



lim Z = 



U'J 1 U'J 1 



Ai^ — 

^ mil V 



2U' Y 2C/'2 r2 



*o = Aie 



iKx 



Ale 



-iKx 



The equations above should be interpreted the following way. Equation (5.2) expresses 
the continuity of the normal velocity at the interface: It imposes the constraint that the 
phase and amplitude of the perturbation of the surface wave (given by A2) is the same as 
the perturbation of the wind (given by ViAi -\-'D2J ) in cluding the component that comes 
from the phase change at the interface. Equation (5J^) is Newton's law, or alternatively 
can be viewed as a statement about the continuity of pressure at the interface; it gives the 
growth of the amplitude A2 due to the out of phase component of the pressure. Equation 
( ^.4[ ) gives the evolution of the vorticity inside the critical layer that determines the phase 
change, and involves the non-linear term "^q xZ^y- The coefHcients Ci, C2, ^^i, ^^2^ € 5ft are 
obtained from the linear theory, and are given in equations ( 4.22 ), ( 4.24 ), ( 4.48 ) and 
(4.49). Ci and I?2 are terms that involve small correction to the real part of C due to 



the gradient of the velocity at the interface U^y and due to the part of the pressure 
perturbation that is in phase with the traveling wave. C2 and Pi involve the part of the 
pressure perturbation that is out of phase with the traveling wave; their sign determines 
the nature of the instability. 



Dropping the non-linear term in eq. (5.4) we obtain J 
fe Reid 1981). Assuming an exponential growth rate Ai - 
be written as: 



iTT^i in the linear case(Drazin 
the above equation can 



e^^a. 



7rC2ai +(7 + «Ci)a2 = , (5.5) 
(Pi - i7rP2)ai -02 =0; (5.6) 
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the growth rate then is given by 



For a positive growth rate we therefore need to have Cal^i < 0. 

We note as an aside that there are several conservation laws at work hcre.| First, the 
vorticity is conserved inside the critical layer {Z)^t — 0, which implies (Z) = since the 
initial Z had infinitesimal amplitude. Second, by noting that (^'o.x^) = iK( J*Ai — JAl), 
one can show that the following laws hold: 

A.{K\A2\'-C2V^(YZ)}^0; (5.8) 

^ IkuM' + \C2V,{Z^)\ ^ -C2V,i.{Z\) ; (5.9) 



|l?iC2((*o + ^(7^r')Z) +P2C2I - Ci|A2p| = yU'C2{Z) . (5.10) 



r„ . 1. 

Recalling that the velocity inside the critical layer is given by \u — C, w] = [^U'^Y + 
e^(l/2[/^'r2 + + • ■ • ,e^*o,2; + ■ • ■] we can identify the first relation (eg to 

correspond to the conservation of momentum. Combining equations (|]^) and ( |5.9| ), we 
obtain the conservation of enstrophy inside the critical layer, 

A.{{U:,'Y + Z)^) = -,.{Z'y). (5.11) 



The third equation ( 5.10 ) can be regarded as a statement of the conservation of energy. 



5.2. Quasi-steady state 

An interesting limit in the set of our equations is when the rescaled viscosity is large 
enough to play a dominant role inside the critical layer. With this assumption we can 



drop the time derivative term in equation (4.13); the functional form of Z then just 
depends on the value of Ai, which then makes J just a function of Ai. More specifically 
we have 

lyZ^YY + *o,x^,y - U'.YZ,, = -f/^'^o.. ; (5.12) 
by letting Ai = i?(r)e*®('^) and £,= Kx + e,we have 

lyZ^YY - 2RKsm{^)Z^Y - KKYZ,^ = 2U'jRKsm{0 • (5.13) 

Let Y = V2'n{R/U'J^^^, Z = y/R/2U^U'jZ{C, rj) and A v^J{K{2Rfl'^)- we then 
obtain 

AZ - sin(0^,^ - r;Z 4 - 2 sin(C) ; (5.14) 



Equation (4.14) then becomes 



^ / e-'^^ZdxdY ^ —R-e:^ / e-'^Zd^dr] 

IK ^~ j-00 J-TT+e 



f In the equations (5.^), ( p.9[ ), and ( p.l[l| ), the integration over x was assumed to be taken 
first. 
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27r 



R-e 



+00 /'+7r 



— 00 — TT 



COS £,Zd£,dri 



+ 00 Z' + TT 



— 00 >.' — TT 



sin C^Zd^drj 



-i — Ai 



00 Z' + TT 



or 



with 



J = -iAi$i(A) 



*i(A) 



27r 



7r /'+00 



sin^Z(i?7c?^. 



$i(A) was first studied numerically by Haberman (1972 ). Its values range from — tt for 
A ^ CXI to for A — > 0. Its asymptotics for A — > and A — > cxd are given by 



$(A) = 



-aiA + a2A3/2 + C'(A2) A<1 



(5.15) 



-7r + feiA-4/3 + o(A-8/3) A>1 
where ai — 5.5151 . . ., 02 = 4.2876 . . . and 61 = 1.6057 .... The derivation for the above 



asymptotics can be found in (Churilov & Shukman 1996). We illustrate the behavior of 
$(A) in figure. (P, constructed by combining numerical evaluation with the aforemen- 
tioned asymptotics. 



cj) -nil 




Figure 1. The Haberman function 
The amplitude equation then can be written as 

A2.T + iCiA2 = C2A1 



[Ki2\A,\)V^ ) 



A2 ^ ViAi + iV2Ai<^> 
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, with 

K{2\A,\)V^J ■ 

We have thus ended up with an ordinary differential equation for the ampUtude of the 
wave. By setting Ai — R{T)exp[iQ{T) — iCiT], and after some algebra, one can show 
that 

(vl+Vl^^^ R^T ^ViC2R^^ , (5.16) 



x>|$$j e,T = C2r'2**, (5.17) 

where $ = (i?<I>).7j. For large R at late times we can conclude that 

\A2\ ^ \Ai\ = R ^ i^T^/^ (5.18) 

and 

e-i/^T-^ (5.19) 

This is one of the basic results of our calculations. As we are going to show later, 
this result becomes valid for later times even for cases for which the rescaled viscosity 
is smaller than one. An important implication of this result is that the amplitude grows 
with an algebraic power instead of the initial exponential variation. Another important 
feature is that the growth rate linearly depends on the viscosity, unlike the case in linear 
theory. (In linear theory, a weak viscosity was giving the same phase change "— itt", but 
the resulting growth rate was independent of i/.) Finally we note that the phase of the 



waves 8, goes assymptoticaly to zero at late times according to equation ( 5.19 ) 

5.3. Numerical Results 
Next we i nvestiga te the weakly non-linear evolution of the wave by solving the set of 



equations ( 15.1^5.41 ) numerically. To solve the advection equation ( |5.4| ), we used a code 
which is spectral in x and finite difference in Y. The domain we used was (—50, 50) in Y 
and (0, 27r) in x. Up to 1024 grid points were used in the Y direction and 63 modes were 
kept in x. The boundary conditions were satisfied to order 1/Y, although the asymptotic 



behavior of Z was taken into account when we evaluated the integral in equation (5.3). 
The code was tested by comparing with a fully pseudo-spectral code as well as with 
already published results. 

Before we begin solving the amplitude equations, however, we re-scale our system so 
that we are left with a minimum number of free parameters. As in the previous section, 
by letting Z ^ U"/VU'Z, Y = v/VTP, T = r/iVTrk), ^ = Kx a-nd u = ly'/VTF we 
obtain the following equation to be solved: 

Z^r + - ^'0,?^,r, - l^'Z,r„ = vE'o.J , 

with 

J^TT / e-'^Zdr^d^. 

Furthermore, by rescaling Ai to |I?ip/|C2p^i, and A2 to |'Di|/|C2p^2, we can always 
scale our system so that 2?i = 1 and C2 = —1- Finally, the coefficient Ci can always be set 
to zero by performing a Galilean transformation (r ^ t + Ci£,) and shifting the critical 
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layer hyY^Y — Ci. The last transformation corrects the position of the critical layer to 
order e, and justifies the neglect of the corresponding terms in §4.2. We are left therefore 
with one independent parameter, 7)2, to investigate; this parameter is a measure of the 
feedback of the gravity wave to the critical layer (I?2 = gives the "free" evolution of a 
critical layer uncoupled from gravity waves). In the following, we examine five cases, for 
five values of I?2 ■ 






(b) 




(c) 


litude 




i change 




Am pi 


/ 
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/ 

/ 
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Time 



Figure 2. The evolution of the wave amplitude \Ai\ (panels [a] and [b]), and of the phase change 
across the critical layer (panel [c]), as a function of time. The dashed line panel [a] gives the 
linear prediction; the solid line in panel [b] is intended to illustrate T^^^ scaling. 







Figure 3. Time sequence of the evolution of the vorticity inside the critical layer for ©2 ~ 0, 
using a grey-scale representation for the vorticity. 



We first examine I?2 = case, e.g. the evolution of a 'free critical layer'. In figure (Qa) we 
plot the amplitude as a function of time; the dashed line gives the linear prediction. The 
amplitude grows exponentially for early times, but as non-linearities become important 
the growth rate decreases, and an oscillatory behavior begins, figure (^jb) shows the same 
thing, but plotted on a log-log scale, and compared with the T"^^^ scaling; we note that the 
agreement is good for later times, although the persistence of the oscillations indicate 
that viscosity is not yet fully dominant. Figure ^p) gives us more information about 
the system; here we plot the phase change 3{— J/^i} as a function of time. The phase 
change (which is responsible for the instability) remains constant, and equal to — tt, in 
the linear regime; in the non-linear regime it drops to zero and then oscillates around 
this value. In figure (^ we display the evolution of the total vorticity Z + r/ inside the 
critical layer, using a grey-scale representation; time frame were taken every 0.5 time 
units (the first panel is at T = 1). The vorticity of the mode rises, as predicted by the 
linear theory, by extracting vorticity from the mean flow. (From eqs. ( ^.9[ ) and (5.11) we 
know that the total vorticity and enstrophy are conserved; therefore the increase of 
must lead to a decrease of ZY.) As the non-linear term becomes important, the mean 
flow advects the vorticity in the vertical direction, leading to the 'winding up' motion 
seen in the panels with T > 2, forming the cat's eyes pattern. The imaginary part of 
the integral J (eq. (|5.3|)) expresses the excess of vorticity of the right side of each panel, 
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from the left. After the first turnover, the vorticity is redistributed inside the critical 
layer more evenly, leading to smaller values of J. Comparing figure (|^) with figure 
it can be seen that deviations from linearity and the —ii: phase change start when the 
vertical motion begins, and the phase change has dropped to zero when the first full turn 
over has been completed. At later times and after a few turnovers, small scale structure 
has been generated in the form of a 'twisted' vortex sheet with sharp boundaries. As 
this procedure carries on, at some point in time, viscosity (now matter how small) will 
become important, and the advection terms in eq. (5.4) will be balanced by the diffusion 
of vorticity; at that point, the asymptotic behavior predicted by the quasi-steady state 
for large times can been shown to hold. 



12 3 





Figure 4. As in figure m, but for V2 = 0.1 



I I -^-i^ 

^^^^^^^^^^^^1 -40 ^^^^^^^^^^^^H -40 ^^^^^^^^^^^^1 



17 3 4 5 6 




Figure 5. As in figure (Q), but for ©2 — —0.1 

Next we investigate the effect of non-zero I?2- The cases we will examine will be for 
2?2 = ±0.1 and ±0.3. We start with the two cases 2?2 = ±0.1; here we note that for 
2?2 = —0.1, linear theory predicts a right traveling wave, but for 2?2 = 0-1 predicts 
a left traveling wave. In figures (^ and (|^) we show grey-scale representations of the 
evolution of vorticity for these two cases. We see that the result is two traveling vortices, 
as predicted by the linear theory, but the phase velocity is decreasing with time and the 
vortices seem to stop at a distance ±7r/2 from their initial position. This is in agreement 
with the quasi-steady limit, which predicts that the phase 6(T') in equation (5.19) should 
decrease as T^^. Small scale differences appear in the vorticity field when compared to 
figure (H), but these seem to be minor. 

By considering lightly larger values of I?2 (= ±0.3), we can obtain an indication of 
how this affects the amplitude of the wave. In figures. (||)-(0) we plot again grey-scale 
representations of the vorticity for four different times for the I?2 = 0.3 and V2 = —0.3 
cases, respectively. The linear growth rate is 1.8 times smaller, so the plots are at later 
times than the previous cases shown. Differences in the evolution of the vorticity can be 
seen especially at the 'saddle point' (the point where the two vortices meet). 

A more quantitative appreciation of the time evolution of the wave amplitude can be 
obtained by directly plotting the evolution of the amplitudes Ai and A2 for the four 
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Figure 6. As in figure (||), but for O2 = 0.3 
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Figure 7. As in figure (|), but for V2 = -0.3 




Figure 8. Evolution of the wave amplitude |y4i| for the four cases, 1)2 ~ 0.1 (panel [a]), —0.1 
(panel [b]), 0.3 (panel [c]), and —0.3 (panel [d]). In each panel, the solid lines show the magnitude 
of the amplitude, while the two sets of dashed lines show the evolution of the real and imaginary 
parts of the amplitude. 
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Figure 9. As in figure (@), but for wave amplitude |y42|. 



cases of non-zero 1^2 we have examined, as shown in figures (§) and (^. In all panels, we 
show both the (unsigned) wave amplitude magnitude, as well as the real and imaginary 
parts of the amplitude; each panel of each figure corresponds to a different value of I?2, 
as indicated. 

Consider first the cases I?2 — ±0.1. The fact that a positive imaginary part is dominant 
in panels [a] of figures (^)-(^), and a negative imaginary part in panels [b] of figures (||)- 
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just expresses the fact that in this case, the two vortices moved by a quarter of 
the wave length and then stopped. We note that the magnitude of the amphtude at the 
point where the oscillations start (e.g., the peak of the 'first bump') does not seem to 
have been affected by the value of I?2- Now consider 2?2 = ±0.3. We mentioned just above 
the appearance of small-scale differences in the grey-scale representations of the vorticity 
field as the magnitude of I?2 is increased; however, the quantitative plots of the wave 
amplitude show that these differences in the small-scale structures do not seem to affect 
the evolution of the amplitudes Ai and A2: for example, in panels [c] of figure (^, we 
see that the value of Ai at which the first bounce of the amplitude occurs does not seem 
to depend strongly on the value of I?2- We conclude therefore that the basic features of 
the weakly non-linear evolution of the unstable mode are not affected by the exact value 
of I?2, at least in the examined range —0.3 < I?2 < 0.3. 



6. Mixing in the cat's eye 

In the previous sections we have examined the cat's eye vortices in the nonlinear 
evolution of the critical layer in the wind coupled to a surface wave. We have found 
that the general flow structure to be similar to those in cases without coupling to the 
gravity waves. Such vortical structures are typical of the weakly nonlinear evolutions of 



parallel flows ( Balmforth fc Piccolo 2001 ). However, their mixing properties may be very 
different (depending on the details of the underlying flow) in spite of great similarity in 
the general flow features. We thus investigate the mixing due to these cat's eyes from the 
weakly nonlinear amplitude equations for a few different cases. 



In our amplitude equation 5.4, the vorticity Z is advected by velocity (it, v) = {Y, — ^oa:)j 
with 'So = Aie^^^ +C.C.. Thus the particle trajectories {x{xo;t),y{yo; t)) satisfy equations 

'i-^-y. (6-1) 

^ = V = aa(t)sm{kx + (t)(t)), (6.2) 
at 

where Xq = (xo,?;o) is the initial particle position and ao{t) and (j){t) are such that 
ikAi = ao(Oe**(*^. 

If Ai is constant in time, equations ( |6.l| )-( |6^ ) reduce to the equations of motion for 
a pendulum, only in this case the particle trajectory is not restricted to a cylindrical 
surface, and can move from eddy to eddy. We can calculate the strain rate of such a 



flow by first linearizing equations (6.1)-( |6.2| ) for an infinitesimal separation between two 
particles Sx, 

^ = (<5x.V)v= f , ^ ■ (6-3) 



Ao, which is the negative determinant of the matrix in equation ( |6.3| ), is interpreted 
as the combination of strain and rotation: Ao(xo) = kaocos{kx{xo;t) + (/>). After time- 
averaging over the tracer trajectories, (Aq) is only a function of the initial position of 
the tracer, and positive (Aq) implies the possibility for a positive Lyapunov exponent for 
that initial position; negative (Aq) implies that rotation is dominant over strain. We have 
calculated the time-averaged (Aq) for each initial position xq in the vorticy for a simple 
case in Balmforth & Piccolo (2001) (figure |l0|(a)). The corresponding plot for the case 



of coupling the shear fiow to the gravity surface waves is shown in figure |To|(b). 
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Fig ure 10. (a): (Xn) for the critical layer without couphng to the gravity wave {T>2 — 
in (Balmforth & Piccolo 2001)). (b): (Ao) for the critical layer coupled to the gravity wave 



(©2 = 0.2) 



The relation between (Aq) and the finite time Lyapunov exponent may not be straight- 
forward, and can depend sensitively on the prescribed flow. In general, some correction 



to (Aq) can be made so th at it is closer to the Lagrangian description ( BofFetta, Lacorata, 
Redaelh fc Vulpiani 200ll). 



A = (Ao) + (Ai), 



Ai = 



xy 



dlpxx dip. 



vv 



dt 



dt dt 



(6.4) 
(6.5) 



where ij) is the stream function of the flow. In our case Ai = 0, and thus we expect (Ao) to 
be a good indicator of the finite time Lyapunov exponent, which we have also calculated 
and shown in figure (p|) for the two cases in figures (p^)a,b. 

From figures (p^)-(pl|) we conclude that the coupling of shear flow to the surface wave 
induces nonlinear evolution in such a fashion that more chaotic mixing may occur in the 
latter case; this conclusion follows from the observation that there appear to be more 
stripes of positive Lyapunov exponent in the driven surface wave case. 

In situations where the passive tracers are weakly diffusive, the asymptotic mixing 
property is determined by the combination of slow diffusion and fast advection. If we 
define 9 as the tracer concentration, we can write down the equation for the weakly 
diffusive tracer in the above flow field, (equations [67 



+ (ao(t) sin(fca:: 



and |6.2D 

-mdy 



1 



= — 
Pe 



(6.6) 



where Pe = UqI/k is the Peclet number, with Uq and / the characteristic velocity and 
length deflned in previous section, and k the molecular diffusivity of the tracer. The 
particle dispersion (variance), often indicative of mixing property of parallel flow, is 
deflned as 

^ {x') - {xf, (6.7) 
where (x) and (a;^) are, respectively, the flrst and second longitudinal moments of the 




concentration field 9 



If the flow is weak and diffusion is strongly dominant over advection, particles undergo 
random walks and their dispersion (cr) increases linearly with the square root of time: 
(7 = (2i/Pe)^''^. On the other hand, if the shear flow is strong and irregular in time, 
the particles will be in a super-diffusive regime during which the dispersion grows faster 
than that for ballistic transport [a ^ t). In cases where the flow is bounded and time- 
independent, the super-diffusive regime eventually gives way to yet another diffusive 
(Taylor) regime with a larger effective diffusivity than molecular diffusion. 

For a time-dependent velocity field, the super-diffusive regime will be the asymptotic 
limit of particle transport, and chaotic mixing may be found instead. To examine how the 
time dependence of the flow affects the particle mixing in our case, we have integrated 



equation (p.q) using a particle method ((Latini & Bernoff 2001) and references therein) 



with ao(t) and (j){t) resembling those from solving the amplitude equations. We place 
10^ particles at an initial position close to the separatrix around the core, solve for their 
positions according to equation (p^) for Pe = 10^, and record their positions from which 



we can calculate the particle dispersion. In our case (equations 6.1 and 6.2), ballistic 
transport is expected (and confirmed numerically) if there is no time variation in Ai. 

Figures ( p^ show the dispersion a versus time for three time dependences of the 
amplitude qq. The solid line is for an exponential growing amplitude which saturates to 
oo '--^ 1 at t '--^ 100, and then oscillates periodically with a fluctuation amplitude of 0.2 and 
a period of 6. The long dashed line is the same except the oscillatory period is 65. The 
dash-dotted line is the same as the first two cases except the oscillation is replaced by an 
algebraic growth proportional to t^/^. We clearly see that the early time dependence from 
simulations is verified as the diffusive regime, where a = (2i/Pe)^/^. We also observe that 
the long time asymptotics for all three cases are at least super-diffusive, with a ~ t^^"^ 
for cases 1 and 2, while for case 3, where the amplitude grows as t'^^^, the dispersion is 
close to cr ~ . 

These results are consistent with anomalous diffusion found in other two-dimensional 
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Figure 12. Particle dispersion for three time- varying amplitudes. 



flows ( Venkataramani, Antonsen & Ott 1998). The temporal periodicity in the flow cre- 
ates the possibility for KAM regions to coexist with chaotic regions, and thus the anoma- 
lous particle transport. However, for the third case, the time variation of the flow is not 
periodic, and there does not exist any KAM regions, and thus the particle transport is 
more super-diffusive. 



7. Discussion and Conclusion 

In this paper, we have extended our earlier linear analysis of the resonant interaction 
between a wind and interface waves to the weakly nonlinear regime. 

Our principal result is the demonstration that the exponential growth of unstable 
resonant waves characteristic of the linear regime transitions to algebraic growth in the 
weakly nonlinear regime, with the wave amplitude growing as t^/^, similar to the cases 
without coupling with the gravity waves. We also flnd that the algebraic growth is linearly 
proportional to the weak viscosity inside the critical layer. 

Before continuing, we return to examine our assumptions which allowed us to carry 
out the weakly nonlinear analysis. As noted earlier, our analysis is carried out around 
the point where the most unstable mode becomes marginally stable (which corresponds 
to, for example, the condition TG = 1/4 for r = 0). In terms of dimensional parameters, 
this means that we conduct our analysis for values of the maximum wind velocity U in 
the range 

Ui<U<U2, (7.1) 
where C/i is found from the criterion that all the perturbations are neutral in the presence 
of surface tension for r = 0: 

and U2 is defined as the upper bound on the maximum wind speed for which only one of 
the unstable modes has crossed the stability boundary. This latter condition arises from 
the assumption of a periodic domain, which leads to a discrete spectrum, and is crucial 
in numerical simulations. 

First we note that equation |7.2| gives a lower bound on maximum wind velocity (C/2) 
that is independent of the length scale (/) of the wind profile. Equation (7.2) thus gives 
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a general lower bound on the maximum wind velocity U for which our analysis applies. 
For the air-on-watcr case, we substitute p2 ~ 1 gmcm"'^, g = 980 cms^^, and a = 
72.8 dynes cm^^ into the above equation, and obtain a lower bound for the wind velocity: 



Ui ^ 0.2 ms ^. For U below this value in air-on- water system, system is stable (Alexakis, 



Young & Rosner 2002). The upper bound on the maximum wind velocity so that our 
analysis applies can be explicitly expressed in terms of size of the periodic domain Ij,, ct, 
g and p2 



U2W ^ max \ ; , \ ; , (7.3) 

' 27rp2^6 V '^T^P'2h j 

where h is of the same order as the wind characteristic length: lb ^ 0{l). For the air- 
water case, we find that the upper bound U2 to be a few times larger than Ui, and thus 
we expect that our results will be a good approximation for winds with velocity up to 
U2 ~ 1 ms""'^. Winds of this magnitude cause ripples on water-surfaces but not any wave 
breaking. 

Within the range of validity of our analysis, we can further examine when the evolu- 
tion transitions from linear to nonlinear regimes, e.g., the transition from exponential to 
algebraic growth. For small density ratio (the air-over- water case), and for high Reynolds 
numbers (for which the results from §5.3 hold), this transition happens when the am- 
plitude is A2 ~ 60. Re-constructing the dimensional amplitude, we conclude that the 
perturbation grows exponentially until its amplitude reaches h/l ~ 60e^ ~ 6 • 10"^, 
where I is the length scale of the wind (see §2; / is of the same order as the most unstable 
wavelength) and e = r = 10"'^ is the air-over-water density ratio. If, on the other hand, 
the rescaled Reynolds number is small enough so that the quasi-steady state approxima- 
tion (§5.2) is appropriate, then the transition will happen when the Haberman parameter 
A is close to one (A ~ 1; see figure 0). This impHes that the transition ampHtude scales 
as h/l ~ 1/ Re, which is still very small since we are in the almost inviscid limit. The key 
point, however, is that in both cases, the amplitude at which the transition in behav- 
ior occurs is extremely small, and most likely is too small to be observed in numerical 
simulations. Thus, any linear growth observed in simulations can only result from wind 
with maximum velocity much greater than 0.2 ms~^. Our results also indicate that the 



energetic estimate provided by Miles (1999) should only work only when the wind is 
much stronger than 0.2 ms~^. 

Furthermore, we note that the linear regime (e.g., exponential growth) would be ob- 
servable if the wind velocity is larger than order one. For example, if 1 ^ ^ r, then 
the linear growth rate will be larger than order r; this implies a different scaling than 
e = r. (The derivation of the new scaling is straightforward, but depends on the how 
the free parameter G scales with e (e.g., G ^ e"); since this is essentially unknown, we 
have chosen not to pursue this calculation further.) In cases for large U , we would then 
expect both the exponential and the algebraic growth rates to be observable in carefully 
conducted laboratory or numerical experiments. 

For astrophysical cases, for which the density ratio is not so small (r ~ 0.1 ~ 0.5), it is 
much easier to capture both scalings (exponential and algebraic) , in either numerical and 
experimental work. That is, following the previous arguments, the transition amplitude 
between the two scaling regimes is h/l ~ O(O.l), and can be easily varied by changing 
the magnitude of the wind. Numerical simulations are therefore obtainable to follow 
the fully non-linear development of the waves; and this is what we intend to do in our 
future work. Furthermore, it remains an open question as to whether the above discussion 
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applies to situations in which the physical domain is not periodic, such as in laboratory 
experiments. 

Finally, we note that we have obtained an interesting result, namely the enhanced 
mixing at the 'saddle points' separating successive cat's eyes within the unstable interface; 
this enhanced mixing (which we studied by means of inserting Lagrangian tracers, and 
observing their evolution) is a consequence of the dynamics near the 'saddle', where 
mixing between two adjacent vortices appears to take place. These islands around the 
separatrices are commonly found in non-integrable Hamiltonian systems, and here they 
serve as implication of chaotic mixing. Results from the particle analysis further confirm 
that chaotic mixing is a consequence of the temporal behavior of the amplitude associated 
with the global background flow. Even though we have assumed much simpler temporal 
behavior for the amplitude in our simulations, the super-difi\ision found in the simpler 
cases affirms that chaotic mixing should be expected as a result of the instability of the 
shear flow coupled to the gravity surface waves. This may imply that the entrainment 
rate of humidity into air could be enhanced by the coupling of weak wind with suface 
waves. 

The calculations presented here will serve as a useful constraint for future numerical 
simulations of wind-driven interface instabilities. Thus, paradoxically, while tests of nu- 
merics involving the linear behavior of such systems can be done, they are challenging 
to apply; in contrast, it should be relatively easy to test the numerics in the weakly 
nonlinear limit. 
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monuclear Flashes at the University of Chicago. We thank N.J. Balmforth, N.R. Lebovitz 
and S.C. Venkataramani for helpful conversations. YNY acknowledges support from 
NASA and Northwestern University, and computation support from Argonne National 
Labs. 
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Appendix 

A. Large G behavior 

We are interested in cases for which the factor G is large. As ah-eady discussed, such 
cases not only correspond to the astrophysical limit of strong surface stratification, but 
also correspond to cases for which the growth rate of (linearly) unstable modes is small, 
i.e., to cases for which our analysis is actually appropriate. 

In order to proceed, we need to adopt a specific wind profile; in what follows, we will 
use a profile of the form U = Umaxi^ — e~^). Ho weve r, we note that our basic results 
hold for more general wind profiles. From equation 3.6 we know that the unstable modes 
will have to be of the same order as G, and therefore K ^ 1; this allows us to write a 



WKBJ expansion for the solution of the perturbation stream function of equation (3.1) 
for the wind. The equation we have to solve for large K is therefore 



by rescaling y Ky, we have 



,vv 



U-C 

1 U^y/K,y/K 



= 0; (Al) 



K"^ U-C 



= (A2) 



or 



-F(y/i^)0 = O, (A3) 

where _F is a slowly varying function of its argument, as illustrated in figure Al. The 
boundary condition at the interface is: 



KC"^ -r[KC'^(t)y\o-CU'\o\ -G(l-r) =0. 



(A4) 



From figure Al it is obvious that the WKBJ approximation will break down at two 




points: The first one is at y = yc^ where the critical layer is located (and where F 
becomes singular); the second one is at y = j/a, where F = 0. For this reason, we will 



Wind-driven gravity waves 

have to decompose the y— axis into three regions: I. < y < yc] II. y 
ya < y- The first-order solutions of the WKBJ equations therefore are: 
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< y <ya; HI. 



(I) 
(11) 

(III) 



B 1 e+Jo"'"'^^', 



A2^sin (^Jlwdy' ~7r/4 



Bo 



—7= 



[Jlwdy'-n/i 



(A5) 
(A6) 

(A7) 



where 



w{y) = 



1 



1 c/. 



.yjK.ylKi 



1 



1 



-v!K 



K'^ 1 - C - e-a/^ 



and the — 7r/4 factor appearing in the solution for Region II is inserted for convenience, 
to be exploited shortly. The coefficients Ai, Bi, A2, B2, are connected through the 
solutions at the points where the WKBJ approximation breaks down, and can be obtained 
using matched asymptotics. Thus, close to y — ya it is well-known that the solution is 
an Airy function; the consequent connection formulae are given below (Olver 1997). For 

y < Va, 

rv 



2A, 



sm 



'w 
243 

'W 



wdy' + tt/4: 



sin [Ii] cos 



wdy' — 7r/4 



— cos [Ii] sin 



wdy' — tt/4 



and therefore 



B2 - 2A3 sin [h] , A2 = -2A3 cos [h] , (A8) 
where Ii — J^" wdy' . The solutions near the critical point y = yc satisfy the equation 

m' 



The solutions of the above equations are given in terms of z — —{y — yc)U'^ /U'^K > by 

/i(z) = .fzJi (2Vi) , f2{z) = .fzNi (2Vi) , (AlO) 

where Ji, Ni are, respectively, the first Bessel and Neumann (Bessel of the second kind) 
functions. The first terms in the asymptotic expansion for y +00 are 



A(z) ^ sin (2zV2 _ ^/A ^ ^^(,) ^ _^ 

Matching with the outer solution, we have 

(j) = V¥A2/i(z) - ^/¥B2/2(^) . 

The asymptotic expansion for z ^ 0^ is 

fi^z + ... 
7r/2 ~ -1 -Fzln|z| + .. . - (1 - 27)z-|- 



(2zi/2 - 7r/4) . (All) 



(A12) 
(A13) 



where 7 is the Euler-Masceroni constant. Thus, we can identify /i with the regular 
Frobenius solution and /2 with the singular Frobenius solution (pi,- Here we assumed 
that Ci is much smaller that (this is something that will be justified a posteriori). 
For y < yc the solutions of equation (K9) can be obtained by making the transformation 
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z e~" z which is equivalent to taking the contour below the critical layer. By doing this 
the first Bessel function transforms to the first modified Bessel function that is growing 
exponentially and is real while the second one transforms to a linear combination of an 
exponentially growing and an exponentially decreasing modified Bessel function, and due 
to the presence of the logarithm it is going to have an imaginary part. Their asymptotic 
expansion for 2/ — > oo is 



h 



(AM) 



2^ ' 20F 
with z = {y — yc)U" /U^K > 0. The inner solution for negative large z can be then written 



as 



Pin 



= Vtt {A2fl - 82/2) 



2 1' 

Matching with the outer solution then gives 



(A2 + iB2) e^^ - B2e-^^' 



2v^ 



\A2+iB2)e •'yc + ±>2e 



or 



A2+iB2 1 -PWy' , B2 1 PWj/' 

(b = e^—=e Jo " + — e ^—=eJo 

^ 2 y/w 2 y/w 

where I2 = Jq" wdy' . Gathering all the terms then gives 



and 



'uy . ^Jaullelmg an uiie ueims uiieii gives 

LI = ^3 cos [h] - i sin [7i] ^ e+^^ , Bi = A3 sin [7i] e"^^ , 
sm{h)e-''e+Io'"'^'"' + cos(7i)e+^=e- T - isin(Ji)e+^=e+.C'"''^' 



The values of </> and c/i^j, at zero arc 

)e^^^ — i sin(/i)e 



: therefore 
)e-'^ + cos(/i 



- cos(/i)e+-^^ + isin(/i)e- 



wind profile (1 



(p\a = A3 [sin(/i 

(p.vla = ^3 [sin(/i)e 

where we have kept terms only to first order in K. For the given 
we have: 

2/e = -Kln{l - C), Va = -[Kln{l - C) - Kln{l + 



h^yc + 0(l/if') = -K ln(l -C)+ 0{1/K) . 



(A15) 
(A16) 

(A17) 
(A18) 

(A19) 

(A20) 

(A21) 
(A22) 
-e--), 

(A23) 
(A24) 



sin(/i) 

Normalizing so that (j)\o - 



■ — , cos(/i) ~ 1, e 
1, wc then obtain 



-I1 



_ sin(/i)e--f2 _cos(/i)e+^= + i sin(Ji)e+^2 
sin(/i)e"'f2 _|_ cos(/i)e+-^2 - i sin(7i)e+-''2 ' 



(A25) 
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or 



= -1 



2sin(/i)e-^2 
sin(/i)e~^2 4- cos(/i)e+-f2 — isiii{Ii)e+^-^ 



(A26) 



The second term in equation (A26) is exponentially small when compared to 1 since 



I2 ~ K; neglecting this term when appropriate then allows ip^y to be written as 

= -1 + 2i sin2(/i)e^2/2 ^ (A27) 

where we have kept only the first term in the expansion of the real and imaginary parts. 
Plugging in this value of (f>,y in equation (A4), we obtain, to zeroth order, 



Co 



G 



l + rK 



(A28) 



which corresponds to the gravity wave in the absence of a wind; At is the Atwood number. 
For our purposes, this is as far as we need to go in analyzing the real part of C. 

We next turn to analyzing the imaginary part of C . To obtain the first order in 
Im{C} = Ci <C Co we have: 



1 



2i^(l - r)CoC, - rKC'oil + (t>.y) = 0, C, = - 



2 1 



-Co{l + (t>,y), 



(A29) 



so that 



sin2(/l)e-2^^ 



(A30) 



Im{Ci} 



AG 



4(1 - r) AG^ V K 



5/2 



(K/AtGY 



VK/iAG), 



2AtG 



(A31) 



We note that (1 — l/^/x)'^'^ is a bounded function smaller than one, and therefore Ci 
has a negative exponential dependence on G. We note further that this exponential 
dependence should be independent of the wind profile, and in a more general case — 
for which U{y) is the wind profile and U^^{c) is its inverse — the growth rate will 
be proportional to d ^ exTp[—2KU^^{c)]; this can be re-written as Gi ~ f{c{K))-^^^, 
with /(c) = exp[— 2J7~^(C)/C] a bounded function and C = Co- The interpretation of 
equation ( |A30| ) is straightforward: it simply states that the growth rate is proportional 
to the negative exponential of the height of the critical layer, as measured in units of the 
wavelength. 
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